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Abstract 

'j^. We derive an exact single-body decomposition of the time-dependent 

"^ I Schrodinger equation for N pairwise- interacting fermions. Each fermion 

C^ ' obeys a stochastic time-dependent norm-preserving wave equation. As 

^ ■ a first test of the method we calculate the low energy spectrum of He- 

lium. An extension of the method to bosons is outlined. 
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1 Introduction 

Solution of the Schrodinger equation for pairwise interacting identical fermions 
is a difficult computational problem with applications in many areas of chem- 
istry and physics. The development of accurate and computationally efficient 
schemes for calculating the ground and excited electronic states of molecules 
is a longstanding goal of theoretical chemistry P . Electron dynamics plays 
an important role in molecular electronics |2] and atomic and molecular dy- 
namics in strong time- varying external fields p. The iV-body problem for 
fermions also arises in shell models in nuclear physics [3]. Exact strategies 



for A^-body problems generally have computational costs which scale expo- 
nentially with the number of particles. Here we show that exact solutions 
of the A^-fermion time-dependent Schrodinger equation can be obtained via 
a multi-configuration Hartree-Fock Ansatz in which the single-particle wave- 
functions for each configuration obey norm-conserving stochastic wave equa- 
tions. Since all properties of the A^-fermion problem can be calculated from 
the exact time-evolving wavefunction, and since the computational costs ap- 
pear to scale favorably with the number of electrons, this method could 
provide a useful alternative to other computational strategies such as time- 
dependent density-functional theory jT] and auxiliary-field quantum Monte 
Carlo 0. 

The technique of decomposing high dimensional deterministic equations 
into lower dimensional stochastic wave equations was pioneered by Gisin and 
Percival |4| who were able to reduce deterministic master equations for the 
density matrix into stochastic equations for a wavefunction. More recently 
the same approach was used to reduce the A^-boson Liouville equation into 
one-boson stochastic wave equations [S]. Similar decompositions have been 
obtained for fermions [H" and vibrations [7] . Unfortunately, the norms of the 
single particle stochastic wavefunctions grow exponentially for the boson and 
fermion decompositions [3 Ej • This is the wave equation analog of the "sign 
problem" which plagues path integral Monte-Carlo approaches |HI. The 
decomposition for vibrations was derived using a stochastic generalization 
of the time- dependent McLachlan variational principle j7], and as a conse- 
quence the equations conserve norm. Here we derive a similar norm conserv- 
ing decomposition for fermions. We demonstrate the use of the method by 
computing the low energy spectrum of Helium. Finally, we explain how the 
same approach can be applied to bosons. 

Before outlining the derivation in section|21we summarise the method here 
for readers who may not be interested in details. In section |3] we explicitly 
prove that the method is exact and that the single body wave equations are 
norm conserving. SectionEldiscusses an application of the method to Helium. 
In section ini we explain how the method can be adapted for identical bosons. 

We consider the general A^ identical particle time-independent Hamilto- 
nian 

N N-1 N 

n^ = Y.Hi^ + i: E vihj) (1) 

where H[i) denotes the single body Hamiltonian instantiated for particle 



i. For electrons in molecules H = — ^^V^/2?Tie — Y,k=i^k&'^/\^ ~ ^k\, for 
example, where the sum is over the nuclei of the molecule. The pairwise 
interaction V{i,j) between particles i and j is represented via 

V{z,j) = J2^u,0,{i)0,U) (2) 

s=l 

as a sum of products of dimensionless one-body (Hermitian or anti-Hermitian) 
operators Og- In section |21 we prove that such an expansion is always possi- 
ble. The coefficients hug have units of energy and may be positive or neg- 
ative. This expansion is developed for the Coulomb interaction e^/|rj — rj| 
in section El (see also Appendix A). Extension of the method outined here to 
time-dependent Hamiltonians is straightforward: simply replace H and Os 
by their time- dependent analogues in (0) below. 

A general initial iV-fermion wavefunction can be written as a weighted 
sum of Slater determinants of A^ single particle wavefunctions. For our pur- 
poses the single particle wavefunctions for a given determinant should be cho- 
sen so that they are linearly independent and normalised but non-orthogonal. 
Each Slater determinant can then be evolved independently. For simplicity 
we now confine our attention to one such initial state 

|vI/(O))=/3A|0i(O)) 102(0))... 10^(0)) (3) 

where A is the anti-symmetrisation operator jU] and /3 is a normalisation 
constant. Here the position of a "ket" in the product indicates which electron 
it refers to, i.e., for \4>i)\4>2) electron 1 is in state \(j)i) and electron 2 is in 
state 102) while for |02)|0i) electron 1 is in state |02) and electron 2 is in 
state |0i). This convention allows us to express some equations more simply 
than would otherwise be possible. 

In our method the exact state |^l/(t)) evolved from (jH)) is constructed from 
the solutions \4>j(t)) of time-dependent stochastic wave equations. Specifi- 
cally, the exact A^ fermion wavefunction is expressed in terms of an average 
M[. . .] via 

\^it))=i3M[A\Mt))\Mt))---\Mm (4) 



where the \4>j{t)) obey Ito-type jTU] stochastic equations 

k^j s=l 

- ^ E E ^s{Os)kOs\<i)j) ]dt + J2 v^^^ {Os - {Os)j) \(p,)dWs 

k^j s=l J s=l 

- yy\uM<pM,)-0^^^^^^^W4^\<Pk)dt 

A..Z.I ^l^^^l<?^^/2(iV-l)Re{(0,|0,)}l^'=^ 

(5) 
for j = 1,...,N. Here we use a notation where {F)j = {4>j\F\(j)j) / {(f)j\(f)j) 
for any single-body operator F. For notational simphcity the exphcit time 
dependence of |0j) and the stochastic random variables dWg has not been 
indicated. [Note that in the case of electrons \(j)j) are similar to the spin-orbit 
single particle wavefunctions of Hartree-Fock.] The symbols dWs{t) represent 
independent normally distributed real stochastic differentials with 

M [dWs{t)] =0 and M [dWr{t)dWs{t)] = 6rs dt. (6) 

The second condition imposes statistical independence of the stochastic dif- 
ferentials. 

Imagine a sequence of time steps all of equal length dt such that t = mdt 
for some integer m. At each time step a set of stochastic differentials is 
sampled from the normal distribution 

P{dW{ldt)) = [l/{2TTdt)]P/^ exp{-dW{ldt) ■ dW{ldt)/2dt} 

where dW^ldt) = {dWi{ldt), . . . , dWp{ldt)) is the vector of stochastic differ- 
entials {p is the number of components of dW). Note that / runs from 1 to 
m. The expectation (^ at any time t can thus be represented in the form 

i^w) = /9 n / dm{idt)p {dw{idt)) A\MmMt)) ■ ■ ■ \<pN{t)) 
1=1'' 

and Monte-Carlo sampling of the integrals then yields the stochastic paths 
generated by Eqs. (0). Each time sequence of sampled stochastic differen- 
tials defines one set of stochastic variables Ws{t) (i.e., Wiener process). Each 
realisation of the set of stochastic variables Ws{t) as a function of time thus 
yields one Slater determinant in the average M[. ..]. Since single particle 



norms are conserved each Slater determinant is equally weighted in the av- 
erage, and error in the mean will scale as 1/vT where L is the number of 
realisations. 

The single particle wavefunctions on the right hand side of (0) are in- 
dependent of the stochastic differentials dWs{t) and so averages such as 
M[F(0i(t), . . . , (l)N{t))g{dWi{t), . . . , dWp{t))] can be calculated via the sim- 
plified formula M[F(0i(t), . . . , 0jv(t))]M[c/(Wi(t), . . . , dWp{t))]. This fact is 
implicit in proofs of norm-conservation and exactness outlined in section |3] 

The fact that all matrix elements (e.g. (0j|O<j|0j)) in the stochastic equa- 
tions involve single particle operators, and the sum over index k ^ j for 
each \(pj), show that the computational costs will scale at least quadratically 
with the number of electrons. For implementations similar to that for He, 
discussed in section the number of terms in the two-body expansion in 
principle scales as the square of the number of electrons (in practice many 
huJs may be small or zero which could improve the scaling of the method), 
and hence evaluation of all {(f)j\0s\4>j) for each j requires A^^ operations, 
making the method scale as 0{N^) overall. The precise scaling is obviously 
model dependent but computational costs should be somewhere in the range 
0{N'^) to 0{N^). Most alternative exact approaches have computational 
costs which scale exponentially with the number of electrons. 

The most important properties of the stochastic decomposition (JH) and (0) 
are its exactness and its norm conservation. 

Equations (0) conserve norm in the mean (i.e., M[{(f)j{t)\(j)j{t))] = 1) 
which gives our decomposition distinct numerical advantages over other de- 
compositions in which the mean norm grows exponentially [H]. In addition, 
our method conserves norm exactly for each individual stochastic realisation 
(see section HI). Note that the norm of A\(f)i(t))\(j)2{t)) ■ ■ ■ I'pNit)) is not con- 
served by our method because the single particle states are non-orthogonal. 
This however presents no problem numerically. 

Using the Ito calculus [TU] we also show in section El that 



ci|^(t)) = (3M 



N~l N 



N 



J2A\Mt))...\d<p,{t))...\Mt)) 



+ E H A\Mt))...\dUt))---\dMt))...\Mt)) 
j=i k=j+i 

= -^HNl'^it)) dt 

h 



(7) 



which imphes that the method is exact for all forms of our equations. 

Explicit time-dependence of the iV-fermion wavefunction is of direct in- 
terest in many chemical problems. Energies can be extracted via the Fourier 
transform of the time auto-correlation function (\l'(0)|\l'(t)). In practice, one 
computes the function 

I{E) = —ReJ^ (M;(0)|M/(t))exp(^— j dt ^ (M/(0)|5(E - 7^;v)|^(0)) (8) 

which will have maxima at the true energies when the end point of integration 
T is sufficiently large. The method therefore also provides access to spectral 
information and in fact it is straightforward to generalize (jHl) so that states 
of specific parity can be extracted. Eigenfunctions can also be obtained. 

2 Single body decomposition of pairwise in- 
teraction 

Consider a general two-body interaction \^(1,2). We will now show that it 
can be expanded in products of one-body interactions according to Eq. Q. 
Let \i) with i = 1,2, . . . denote a complete basis of the one-body space. Then 
\h]i2) = 1^1)1^2) for ii,i2 = 1,2,... will be a complete basis of the two-body 
space. Here again we employ the convention that the position of a "ket" in 
a product indentifies the fermion. It follows then that we may represent the 

interaction via 

00 00 00 00 

r(l,2) = E E E E \^l;^2){^u^2\Vil,2)\J^;J2){J^;J2\ (9) 

«1 = 1 12 = 1 jl = i J2 = l 

where we have inserted closure relations for the two-body space on either 
side. 

Define a bijective application a : N^ —>■ N which maps each couple of 
integers {i,j) in a unique integer a{i,j). Then we can introduce new com- 
posite indices Ui = o'{ii,ji) for body 1 and cr2 = cr(^2; J2) ior body 2 with ai 
and (72 taking integer values 1,2,.... We may then define matrix elements 

V.,,., = (2i;22|^(l,2)|ji;j2) 

which are symmetric under the interchange of ai and o"2. This symmetry 
refiects the indistinguishability of the particles. Diagonalising V then gives 

00 

"^'o-i,cr2 — Z_^^^sQai,sQa2,s (10) 



where %uJs are the eigenvalues and Q^^g are the dimensionless matrix elements 
of the orthogonal transformation. With a slight change of notation and using 
the inverse of the mapping ai = cr(zi, ji) we may then write 

Q.us = {^l\Os\jl) (11) 

which defines the one body operator Og- Since 1^(1, 2) is Hermitian it follows 
that each Og must be either Hermitian or anti-Hermitian. The eigenvalues 
huJs may be positive or negative. 

Substituting (HID into (HUD, and (HOI) into (0) gives 

oo oo oo oo oo 

il=l 12=1 Ji = l J2 = l s=l 

oo oo oo oo oo 

= Jl^^^Yl E E E Ki;^2)(«i;«2|0,(i)0,(2)|ji;j2)(ji;J2|. 

S=l ii=li2 = l ji = l J2 = l 

Finally removing the closure relations gives 

oo 
s=l 

which is the desired expansion. 

In practice a finite basis set is more practical than a complete one but 
the same considerations apply except that the sum will terminate at some 
finite value p. 

3 Derivation of stochastic wave equations 

We originally derived the stochastic decomposition discussed above using the 
stochastic McLachlan variational principle developed in Ref. j7]. Here we 
present a more direct argument. For simplicity we initially focus on just two 
fermions with the simplest possible interaction. Consider then the restricted 
two fermion Hamiltonian 

H2 = H{l)+H{2) + hujO{l)0{2) (12) 

and a normalised initial wavef unction of the form 

\^{o)) = f3{\Mom2m-\Momim) 



where (3 = 1/J2{1 — |(0i(O)|02(O))p) is a normalisation factor, and |0i(O)) 
and 102(0)) are normalised but non-orthogonal states, i.e., 

(0i(O)|0i(O)) = (02(0)102(0)) = 1 and (0i(O)|02(O)) ^ 0. 

Note the antisymmetric form of the initial wavefunction. 

We wish to find stochastic equations for |0i(t)) and |02(i)) such that the 
exact solution |\l/(t)) of the Schrodinger equation 

d\-^{t)) = ~n2\-^{t))dt 

can be written as the expectation value 

|vI/(t))=/?M[|$(t))] (13) 

of the antisymmetric stochastic vector 

|$(t)) = |0l(t))|02(t))-|02(t))|0l(t)). (14) 

We will also require that the stochastic wave equations conserve norm 

(0l(t)|0l(t)) = (02(t)|02(t)) = l. 

To achieve norm conservation the single fermion wavefunctions must sat- 
isfy the condition 

d ((0i|0i)) = (t/0.|0i) + (0iM0i) + (#i|#i) = (15) 

for i = 1,2. Since d\(f)i{t)) will have a term proportional to a change dW{t) 
in a stochastic process W{t) with M[dW{tf] = dt and M[dW{t)] = , 
there will naturally be terms proportional to dW{tY (which is of order dt) 
in condition (jlSp . Hence it may prove useful to have a term proportional 
to dW{ty in d\(j)i{t)) in order to conserve norm. Our wave equations should 
therefore be of the form 

rf|0i) = \vi)dt + \ui)dW + \wi)dW'^ (16) 

where all quantities depend on the time t. With this form of the stochastic 
differential (i|0j), condition (fT3j) can be written as 

2Re {{(pi\vi)} dt+2Re{{(pi\ui)} dW+{2Re{{(l)i\wi)} + {ui\ui)) dW^ = (17) 



and the coefficients of dt, dW and dW^ must independently vanish. 

In order to reproduce the interaction term of Hamiltonian (jl2|) we must 
have a term in \ui) which is proportional to O|0j). To make the coefficient 
of dW vanish in Eq. (fTTj) it would thus be sufficient to choose 

|Mi) = v/^zZ7(O-(O),)|0,) (18) 

eliminating one of the unknowns in Eq. (fT^ . Here (0)j = {(f)i\0\(f)i) / {(f)i\(f)i) 
where we keep the factor of {(pi\(t>i) explicit even though it is unity. 
To make the coefficient of dW'^ vanish in Eq. (J17p we can choose 

since the \wi) terms were included precisely for this purpose. Clearly 0i and 
02 must be non-orthogonal initially and a declining overlap will cause an 
increase of (|19p for each mode thereby restoring the overlap. 

Finally, we need to find |f ,). Clearly, there should be a term like —{i/h)H\(j) 
to reproduce the single particle terms of Hamiltonian (fT^ . There could also 
be a term like O|0j). So assume that \vi) will take the form 

1^;,) = -{t/h)H\<p,) + aM + b,0\<p,) (20) 

where a, and 6j are unknowns. To make the coefficient of dt vanish in Eq. (fTTjl 
it is necessary that Re{aj} + Re{6j}(0j|O|0j) = 0. Hence we can probably 
set the real parts of a^ and 6j to zero. To determine their imaginary parts we 
consider the expectation of the differential of the vector (J14p which, because 
of condition (jT^. must be equal to the differential of the vector |^I/(t)), so 
that one has 

d\^{t)) = (3M[\dMt))\Mt)) + \Mt))\dMt)) + \dMt))\dMt)) 

- \dut))\Mt)) - iMmdMt)) - \dMmdMm- (21) 

Replacing the differential terms d\(f)i) in the right-hand side of the previous 
equation with expression (fT^ and making use of the results (fTHj) and (fT^ as 



well as of Ansatz (|2(Jj). after some algebra one obtains 



+ (ai- 



-{i/h){H{l) + i/(2))|$(t))rft - iuO{l)0{2)mt))dW^ 
a2)mt))dt - iuj{0),{0)2mt))dW^ 



+ V^KJiOil) + Oi2)-{0)i-{0)2)Mt))dW 

+ {bidt + iLu{0)2dW^) |O0i(t))|02(t)) 

+ (b2dt + iu;{0)idW^) \Mt))\OMt)) 

- lb,dt + iu{0)2dW^) |02(t))|O0i(t)) 

- {b2dt + zu;{0)idW^) \OMt))\Mt)) ■ 

(22) 
Using condition ^ and the facts that M[dW] = and MldW^] = dt, and 
assigning 

ai = a2 = ^(0)1(0)2, bi = -iuj{0)2 and 62 = -^^{0)1, 



we then find that Eq. ^ reduces to d\^{t)) = -{i/h)n2\^{t))dt which is 
the exact Schrodinger equation in differential form. Hence we have found 
exact stochastic wave equations of the form 



d|(/.i) 



d\^2) 



-'-H\ct>,) - iuj{0)20\ct>,) + |^)(O)i(O)2|0l) ) dt 



+ V^iu{0-{0)i)\(pi)dW -\uj\ 



nwii 



) (OtO)i-l(O), 



2Re{(0i|02)} 



'-H\<P2) - Z^(O)iO|02) + y (O)i(O)2|02) ) dt 



+ ^/^i^{0-{0)2)\(t)2)dW -\uj\ 



'J2 (P2 



) {0^0)2 -\{0), 



2Re{(0i|02)} 



2)dW^ 



\(j)i)dW^ 



(23) 

which conserve norm by construction. Since terms of order dW^ and higher 
are of no importantance and since the average of dW'^ is dt, it is possible to 
make this replacement in Eqs. (j^Hj) with no loss of accuracy or generality JU] 



10 



giving 



) = (^^l.H\^^)-^co{OhO\^,) + '-^){0),{Oh\<Pl) 



dt 



+ V^iu{0-{0)i)\(pi)dW-\uj\ 



i|0i) (OtO)i-|(0)i 



2Re{(0i|02)} 



\4>2)dt 



d\ 



■'-H\ct)2) - iuj{0)rO\ct)2) + y (O)i(O)2|02) ) dt 



+ V^iu{0-{0)2)\(t)2)dW -\uj\ 



nw2 



) (OtO)2-|(0), 



2Re{(0i|02)} 



\<Pi)dt. 

(24) 

Generalisation of ()24|1 to the full pairwise interaction gives a special case 
of (JH) and (jni). We thus proceed directly in the next section to consideration 
of the iV-fermion problem with full pairwise interaction. 



4 Exactness and conservation of one-body norm 

Consider conservation of norm first. To be norm conserving Eq. (0) must 
satisfy the constraint 

rf((0,(t)|0,(t))) = (#,-(t)|0,(t)) + (0,(t)|c?0,(t)) + (rf0,(t)|#,(t)) = 

for j = 1, . . . ,N or equivalently that 

dM[{(j)j{t)\(f)j{t))] = and dM[{(j)j{t)\(j)j{t))^] = 0. 

Substituting (jSJ in dM[{(pj{t)\(t)j{t))] gives 

ME \u,\{<P,\<P,) ((OtO,), - |(0,),f ) {dW^ - dt)] 

s=l 

which vanishes. Similarly, 

dM[(0,(t)|0,(t))2] = M[2(0,(t)|0,(t))rf((0,(t)|0,.(t))) 

+ 2|(0,-(t)|(i0,-(t))p + (0,(t)|#,(t))2 + (#,(t)|0,.(t))^ 

which then gives 



M 



2 E |^.|(0,|0,f ((OlO.), - |(0.),f ) (dW^,2 - dt) + 0{dt') 



s=l 



11 



which vanishes as (it ^ 0. Hence norm is exactly conserved for individual 
stochastic realisations as well as in the mean. 

Now consider the issue of exactness of the decomposition. Substitut- 
ing (jSJ into Eq. (|7j) we see that the term of (0) proportional to \us\ makes no 
contribution because the Slater determinants have two identical single par- 
ticle orbitals and hence vanish. The term of (0) proportional to dWg makes 
no contribution to the first term of ((Tj) because M[(iiy5] = 0. The first three 
terms of (0) contribute 



M 



N-l 



--EA|0i)...|/70,)...|0^) 
N-1 N p 

+^E E J2^s{0s)jmkA\h)--.\M 

j=l k=j+l s=l 
N-1 N p 

-^E E E^s{Os),A\<i,,)...\Os<i>k)...\ 

j=l k=j+l s=l 
N-1 N p 

-^E E E^«(Os)fcA|0i)...|o,0,)... 

j=l k=j+l s=l 



(25) 



W) 

/)7V) 



dt 



to the first term of ((Tj). The non- vanishing contributions of the second term 
in Eq. ((Tj) are 



M 



N-1 N p 

-^E E J2''sA\(Pi)...\o,<p,)...\o,<Pk)...\<pN) 



j=l k=j+l s=l 
N-1 N p 

-^E E Y.^s{Os),{Os)kA\<p,)...\M 

j=l k=j+l s=l 
N-1 N p 

+^E E T.^s{Os)jA\(p,)...\Os(Pk)...\(pN) 

j=l k=j+l s=l 
N-1 N p 

+^E E J2''s{0s)kA\(p,)...\0s(P,)...\(PN) 

j=l k=j+l s=l 



(26) 



dt. 



The first terms of ((^ and (PUJ) combine to give —{i/h)HN\'^(t))dt. The 
second terms of ()25p and ()26|) cancel as do the third and fourth terms of 
the respective equations. These considerations thus show that (i|\l/(t)) = 
— {i/h)7ii\f\'${t))dt and hence that the method is exact. 



12 



5 Application to He 

Here we apply Eqs. (0} and (0) to the problem of calculating the low energy 
spectrum of Helium as a first test of the method. Clearly we need a basis 
in which to represent the single electron wavef unctions. We also need a 
decomposition of the form ((2)) for the Coulomb interaction. 

We choose to represent the single particle wavefunctions in a finite ba- 
sis of states |?7,,/, m,r) = \ipn,i,m) <S> \t) where ipn,i,m are the exact orbital 
eigenfunctions of the He"*" ion and |r) = |±) for r = ± are the spin-1/2 
eigenstates. Here n = 1, 2, . . ., / = 0, ... ra — 1, and m = —I, . . . , 0, . . . / 
are the allowed values of the quantum numbers. That is, the orbital parts of 
these basis functions are exact eigenfunctions of the single-body Hamiltonian 
—h V^/2r7ie — 2e^/r with eigenvalues En = —2jr? in atomic units (i.e. ^ = 1, 
rUf. = 1, and e = 1). The functional forms in the coordinate representation 
are 

Wn,l,m) = Rn,l{r)yi,mi^A) (27) 

where yj,m(6', 0) are the usual spherical harmonics (i.e., eigenstates of angular 
momentum) and the radial functions Rn^i{r) are 



4 

RnA^) = -^ 



n^ 



where L^{x) are associated Laguerre polynomials defined as 

Note that this definition of the associated Laguerre polynomials while consis- 
tent with standard mathematical usage ^T] differs from those used in some 
standard physics texts ^12j . 

In this basis the coefficients of the single particle wave functions are de- 
fined via 



c, 



■n,L,rW = (ri,/,m,r|0i(t)) 
c2L,rW = (ri,Z,m,r|02(t)) 



and the components of the full wavefunction \l/(t) in this basis will be defined 



as 



(~'n,l,m,T;n',l',m.',T'\t) — jjM C^^i^rn,Ti^)'^n',l',m',T'V^) ^n' ,1' ,m' ,T'i^)^n,l,m,Ti^) 

13 



which obviously incorporate the correct antisymmetry. 

The basis (P7|) is also used to expand the two-body interaction e^/|ri — r2| 
in accord with Q. Specifically, we calculated the matrix elements 



{H;jl\Vil,2)\i2;j2) = (^ni,/i,mi(l)V^n2,«2,m2(2)|- 



|ri - r2 



#n;,i;,m'i(l)^^,«^,m^(2)) 



for which we provide formulas in Appendix A. Here ii{ni, li, rrii), ji{n2, h, ""^2), 
^2(^1, 1'l, m[) and j2(^25 ^25 ^2) ^^^ composite integer indices ranging from one 
to infinity (or the maximum number of elements in the basis set). The pro- 
cedure outlined in section |21 was then performed numerically to obtain the 
expansion of the two-body interaction. If due to truncation of the basis set 
1 < HiJiii2iJ2 ^ K, then the number of terms in the decomposition Q is 
p = K^. In practice it is convenient to consider some n^ax from which it 
follows that K = nmax{nmax + l)(2n„aa; + l)/6. 

The stochastic equations for the coefficients thus take the form 



2i 



(1) 



+ -,dt ci:>it) 



n^ 



dc^nUAt) = (J2{0s)l [-{Os)2UJsdt-V^uJ-sdWs{t) 

\s=l J n',l',m' 

P -J- - . 



I0i)[(ola)i-Ka)ip] . (2) 



=1 



2Re{(0i|02)} 



dt\ Ci(t) 



<U(i)= E(0.)2 



p 

— ^ 

s=l 



■{Os)iOJsdt - ^/-iujsdWs{t) 



2i 



(2) 



+ ^dt ci:i(t) 



n^ 



\s=l ) n',l',m' 

(02102) f(OtO,)2-|(0,)2p^ 



E 



UJ^ 



.(1) 



dt C„; (t) 



(28) 



2Re{(0i|02)} 
where we calculate expectations via formulas such as 

(0l|O,|0i)=E E E c3,m,rW("''^'"^l<^^l^''^''"^')cS,_„,^^(t). 
T n,Lmn' ,1' ,m' 



An initial state consisting of random mixtures of Is and 2s He^ basis func- 
tions for each electron was chosen. We chose a basis set with Umax = 4 to 
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perform the calculations. Equations (j28|) were solved using an order 4.5 vari- 
able time-step (i.e., adaptive) Runge-Kutta method which has been specif- 
ically developed to solve such stochastic differential equations [TU [T3] . 




Figure 1: Re (^(0)|^(t)) vs. t 



A detailed discussion of the computational method will be presented else- 
where pn] • In Figs, n and |21 we plot the real and imaginary (respectively) 
parts of {ip{0)\ip{t)) against time in atomic units for an exact propagation of 
the initial state (solid curve) and for the solution obtained via Eqs. (j^ for 
200000 realisations (dashed curve). The agreement is satisfactory although 
the calculation has not completely converged. In Fig. |21 we show the energy 
spectrum calculated via Eq. (jS)) for the exact and stochastic wave solutions. 
Again agreement is good with the stochastic calculation reproducing all en- 
ergy levels.. 



6 Extension to bosons 

Stochastic decompositions for pairwise interactions are also of interest for 
bosons in the context of Bose- Einstein condensation [Sj. Here we show that 
our approach can be adapted to bosons as well as fermions. 
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Figure 2: Im (^(0)|^(t)) vs. t 



3 - 




Figure 3: He energy spectrum 
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To begin with consider the case of two pairwise interacting bosons. Again 
we break the general initial wavefunction into a sum of symmetric states of 
the form |^(0)) = /9(|0i(O))|02(O)) + |(/)2(O))|0i(O))). We need an exact means 
of propagating the single particle states individually such that 

\^it)) = (3M[\MmMt)) + \hmMm m 

where \4>i{t)) and \4>2{t)) satisfy norm-preserving stochastic wave equations. 
To do this we add a fictitious subsystem of two spin-1/2 degrees of freedom 
with null Hamiltonian and anti-symmetric state (l/v^)(|+)|— ) — |— )|+)) to 
our problem. We thus have a total wavefunction 

l^jUt)) = {p/V2)M[\Mt))\Mt)) + IMtWim ® (l+)|-) - |-)l+)) 

= (/3/x/2)M[(|0i+(t))|02-W) - I02-W)l0i+W)) 

- (l01-W)|02+(t))-|02+(t))|01-W))] 

(30) 
where |0j±(t)) = |0j)|±) for i = 1,2. This wavefunction is a sum of two 
antisymmetric states. It is thus clear that solutions of (|3U|) can be obtained 
by determining the time evolution of two-particle antisymmetric states 

101., (t)) 102.. (t)) - |02..(t))|01.,(t)), (31) 

for 0"i,cr2 = ±, which can be obtained with the method for fermions out- 
lined above. Hence we can obtain (jHUj) at the cost of including an extra 2 
component spin to each single particle state. From ()30|) we can get (j^^ by 
projecting out the fictitious part of the solution via 

I^W) = ^((+l(-|-(-|(+l)l^/.c.(t)). 

Hence the 2-boson problem can be solved using the 2-fermion formalism at 
the expense of doubling the number of equations. 

For an arbitrary number of bosons A^ we wish to find a stochastic decom- 
position 

|vI/(t))=/5M[5|0i(t))|02(t))...|0iv(t))] 

for dynamics generated by Hamiltonian (^ where S is the symmetrisation 
operator. Let 

\a) = aAalcTi) ...\aN), 
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where \aj) denotes one of a set of A^ spin states, and Aa is the anti-symmetrisation 
operator on this space. Here a is a norniahzation constant. If the spins again 
have a null Hamiltonian then we may define a fictitious dynamics 

|^/.ct(t)) = |^(t))|a) = a(3M [SAa\<Pi.,) . . . |0^.^)] 

where each of the A^! terms Aa\(j)iai) ■ ■ ■ \(t>Nai^) in the symmetrisation sum is 
antisymmetric. Here |0jo- ) = \<Pj)\<^j) and so by adding an extra fictitious 
spin with A^ allowed states we can convert the problem into fermion form. 
Application of the fermion method is then straightforward and the boson 
wavefunction can be extracted in the end by projecting out the fictitious 
spin state via |\E'(t)) = {a\^ fidit)) . Because of the need to introduce a 
fictitious spin the computational costs of the boson method scale between 
0{N^) and 0{N^) depending on the nature of the interaction. 

7 Summary 

We have shown that the time-dependent quantum A^-body problem for pair- 
wise interacting fermions can be exactly decomposed into A^ one-body prob- 
lems each of which obeys a stochastic norm-conserving wave equation. Our 
approach improves on previous decompositions [HI because the single par- 
ticle equations conserve norm and thus are much more stable numerically. 
Use of the method was demonstated by calculating the low energy spectrum 
of Helium. We have also explained how the approach can be extended to 
bosons. 

The authors acknowledge the financial support of the Natural Sciences 
and Engineering Research Council of Canada. 



8 Appendix A 

Consider the single particle Hamiltonian H = — /;,^V^/2r?7,e — Ze'^/r which 
has Hydrogen-like eigenfunctions of the form ()27|) with 



«"^^«4\-^^f:^--"(f^)'^^--(f^) 
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in atomic units with associated energies En 
shown that 



X 



X 



X 



X 



X 



+ 



(^ni,ii,mi(l)V'n2,«2,m2(2)| 



(2/i + l) (2/^ + 1) 



ri - r2 



—Z'^/2n^. It can then be 



K -h- 1)!K - I'l - 1)!(^2 -k- 1)!K -1'2- 1)! 



16^(2Zi + l)(2/2 + l)\ 

oo i r /' / / 

Y^ V^ (_T\m '■I '' 'l 

l=Om=-l 



im + h)\{n[ + l[y.{n2 + k)\{n', + l',)\ 



ni\ m mi 



I'o 



rrir, 



-m 7712 



l[ I h 





'2 ' '2 




il-«i-l n[~l[-l m-h-l "2-'2~l 



E E E E (-1)^^+'=^+'=^+'=^ 

fci=0 k[=0 k2=0 kL,=0 



, {h + I'l + l2 + l2 + ki + k[ + k2 + k'2 + 4)! 



h\k[\k2\k'2\ 

■n,2 + ^2 \ I n'2 + ^2 



I ni — /i — 1 — A;i J \ 77[ — l[ — 1 — k'l J \ 712 — I2 — ^ — k2 I \ n'2- 
(2/ni)'=i+'i+2(2/n;)^i+'i+2(2/n2)^2+'=^+2(2/n'2)'^2+'2+2 

F(l, /i + /; + /2 + ^2 + ^1 + ^'1 + A;2 + A^a + 5; 



I'o 



Kn 



l + h + l'i + ki + k[ + 3 
l + h + l[ + ki + k[ + 4; 



1/ni + 1/n'i 



1 



/ + /2 + ^2 + ^2 + A;2 + 3 

/ + /2 + ^2 + ^2 + /^2 + 4 



1/ni + l/n'i + 1/^2 + l/n'2 
F(l, /i + /; + /2 + ^2 + ^1 + ^'1 + A;2 + ^2 + 5; 

1/712 + l/n'2 



1/ni + l/n'i + 1/^2 + l/n'2 
where F{a, b; c; (i) is the hypergeometric function |llj. Here 



Jl J2 J 
mi 7712 ^^ 



{ji,J2,mi,7n2\j,7n) 



denote the Clebsch-Gordon coefficients. We have used the convention of 
Ref. ^2] in which 



Jl J2 J 

mi m2 ^« 



5. 



m,mi+m2 



i2j + l)ABY: 

n=0 



(-1)' 



(32) 
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with 

(J+J1+J2 + I)! 

5 = (ji + mi)!(ji -mi)!(J2 + m2)!(J2 -"^2)!(j +"^)Kj -"i)' 

Cn = iji + J2 - j - n)\{ji - mi - n)\{j2 + 1712 - n)\{j - J2 + mi + n)\ 

where it is understood that the sum in ()32|1 truncates when factorials have 
negative arguments. 
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